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Abstract 

We investigate the positional order of the two-dimensional hard disk model with 
short-time dynamics and equilibrium simulations. The melting density and the crit- 
ical exponents z and n are determined. Our results rule out a phase transition 
as predicted by the Kosterlitz-Thouless-Halperin-Nelson- Young theory as well as a 
first-order transition. 
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Usually, numerical measurements of critical exponents are carried out from 
simulations in equilibrium using Monte Carlo (MC) techniques. Such simu- 
lations at or near the phase transition point, except for some special cases 
[1], are affected by the critical slowing down. Recently, Janssen, Schaub and 
Schmittmann [2] proposed an alternative, which allows their determination 
from the short-time dynamics. They discovered that a system with non-con- 
served order parameter (model A [3]) quenched from a high temperature state 
to the critical temperature shows universal short-time behaviour. This sets in 
after a microscopic time scale t mic , during which the non-universal behaviour 
is swept away. Starting from a small value of the order parameter mo the order 
increases with a power law M(t) ~ mo t 9 , where 9 is a new dynamic exponent. 
This short-time behaviour was supported by a number of MC simulations [4,5]. 
These investigations also provide a possibility to determine the conventional 
critical exponents. The exponents can be calculated from the time evolution 
of the second moment of the order parameter, the cumulant and additional 
observables, which evolve also with a power law. Since the simulations are per- 
formed in the early part of the dynamics, this method may eliminate critical 
slowing down. 

While universal short-time behaviour was first seen when starting from un- 
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ordered states, short-time dynamical scaling can also be found starting from 
the ordered state (M(t = 0) = 1). In the latter case, no analytical calcula- 
tions exist, but MC simulations have been performed [6,7,5]. Again, one can 
use the short-time behaviour to calculate the critical exponents, except for 
the new exponent 6. In this letter, we use the dynamic relaxation of the two- 
dimensional hard disk model starting from the ordered state to calculate the 
critical exponents for the positional order. Also this is the first time that the 
dynamic evolution in the early time is studied for a non-lattice model. 

The nature of the two-dimensional melting transition is a long unsolved prob- 
lem [8,9]. The Kosterlitz-Thouless-Halperin-Nelson- Young (KTHNY) theory 
[10] predicts two continuous transitions. The first transition occurs at temper- 
ature T m when the solid (quasi-long-range positional order, long-range orienta- 
tional order) undergoes a dislocation unbinding transition to the hexatic phase 
(short-range positional order, quasi- long-range orientational order). The sec- 
ond transition is the disclination unbinding transition (at 7]) which transforms 
this hexatic phase into an isotropic phase (short-range positional and orienta- 
tional order). There are several other theoretical approaches to the transition. 
An alternative scenario has been proposed by Chui [11]. He presented a the- 
ory via spontaneous generation of grain boundaries, i.e. collective excitations 
of dislocations, and predicted a conventional first-order phase transition from 
the solid to the isotropic phase. In this case, a region exists where both phases 
coexist instead of a hexatic phase. Even for the simple hard disk system no 
consensus about the nature of the transition has been established. 

A number of simulations of the two-dimensional hard disk model in equilibrium 
have been performed. The melting transition of the hard disk system was 
first seen in a computer simulation by Alder and Wainwright [12]. They used 
a system of 870 disks and molecular dynamics methods (constant number 
of particles N, volume V and energy E simulations) and found that this 
system undergoes a first-order phase transition. But the results of such small 
systems are affected by large finite-size effects. Recent simulations used MC 
techniques either with constant volume (NVT ensemble) [13-16] or constant 
pressure (NpT ensemble) [17,18]. The analysis of Zollweg and Chester [13] 
for the pressure gave an upper limit for a first-order phase transition, but 
is compatible with all scenarios. Lee and Strandburg [17] used isobaric MC 
simulations and found a double-peaked structure in the volume distribution. 
Lee-Kosterlitz scaling led them to conclude that the phase transition is of first 
order. However, the data are not in the scaling region, since their largest system 
contained only 400 particles. MC investigations of the bound orientational 
order parameter via finite-size scaling with the block analysis technique of 
16384 particle systems were done by Weber, Marx and Binder [14]. They also 
favoured a first-order phase transition. In contrast to this, Fernandez, Alonso 
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and Stankiewicz [18]Q predicted a one-stage continuous melting transition, i.e. 
a scenario with a single continuous transition. Their conclusions were based on 
the examination of the bond orientational order parameter in very long runs 
of different systems up to 15876 particles and hard-crystalline wall boundary 
conditions. Mitus, Weber and Marx [15] studied the local structure of a system 
with 4096 hard disks. From the linear behaviour of a local order parameter they 
derived bounds for a possible coexistence region. Finally, we showed by MC 
simulations in the NVT ensemble [16] that a one-stage transition can be ruled 
out. The results of the bond orientational order parameter were compatible 
with the KTHNY predictions or a weak first-order transition. 

In this letter, we present results for the positional order of the hard disk model 
obtained through MC simulations near the melting transition p m and in the 
solid phase (p > p m ). First we use traditional simulations in equilibrium (NVT 
ensemble) to locate the melting density p m and to calculate the critical ex- 
ponent rj. After that, we investigate the dynamic relaxation at and above p m 
(model C) and calculate the critical exponents z and rj from the short-time 
behaviour. We use always periodic boundary conditions and a rectangular box 
with ratio v3 : 2, which is necessary since we start the relaxation from the 
ordered state. The disk diameter is set equal to one so that the lengths of the 

system are given by L and 2L/v^3, where L = J y / 3A^/2p. For the simulations 
of the dynamic relaxation we use the conventional (local) Metropolis algo- 
rithm and choose the new positions of the particles (for symmetry reasons) 
with equal probability within a circle centered about its original position. Sim- 
ulations in equilibrium are performed with a non-local Metropolis updating 
algorithm [20]. In this case, the new positions of the particles are chosen as 
usual within a square. 

The fcth moment of the positional order parameter ?/> pos is defined as 




where G is a reciprocal lattice vector and r*j denotes the position of particle 
i. The magnitude of G is given by 2ir/a, where a = \jlj\pip is the average 
lattice spacing. The orientation of G was defined in two different ways, which 
lead to two different Vw- I* 1 a ^ TS ^ definition we fix the direction of G to that 
of a reciprocal lattice vector of the perfect crystal (which are unique due to 
the boundary condition of a rectangular box of ratio 2 : y/3). The reason is 
that large crystal tilting (with the given boundary conditions) is not possible 
if we are in or near the solid phase, while small fluctuations of G are ignored. 



For a critical discussion of this work see Ref. [19]. 
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However, we use a second definition, where we determine the orientation of 
the crystal from the global bond orientation 



i N i Ni 

^ = ^E^E^ P (6i^). (2) 

1 = 1 1 j= 1 

The sum on j is over the Ni neighbours of the particle i and 0y is the angle 
between the bond formed by particles i and j and an arbitrary but fixed 
reference axis. Neighbours are obtained by the Voronoi construction. Since we 
try to determine the global orientation of a single configuration, no average is 
taken in Eq. (2). With ip Q ~ exp(6ia), we define the angle a (0 < a < tt/3) 
and therefore the orientation of the crystal and G. Of course, the definition of 
G is only valid in a strict manner if the bond orientational correlation length 
^6 is infinite (i.e. if we are in the hexatic or solid phase). Simulations show 
that this is probably the case [16]. However, if £g is finite, we have at least the 
case that £g ^ L, since we simulate near the melting transition. Thus we have 
also in this case some kind of global orient at ionT^]. The simulations showed 
that the difference between the two definitions is negligible in the solid phase, 
while it is getting important below the transition point. 

The positional correlation is quasi-long-ranged in the whole solid phase. There- 
fore, the fourth-order cumulant 




(measured in the equilibrium) is independent of the size of the system for 
P > Pm- We use this finite-size scaling (FSS) behaviour to locate the melting 
density p m with simulations in equilibrium. Our results for the dependence of 
U on L for N = 32 2 , 64 2 and 128 2 are shown in Fig. 1. We use the second 
definition of ip pos (varying orientation of G), since most of the simulations 
are performed below p m . However, the other definition leads to comparable 
results. The figure shows that the scale invariance of U yields a melting density 
of p m ~ 0.933 (in reduced units). Since there is a tendency of the slope to 
decrease with larger systems, scale invariance may actually take place at a 
slightly larger density. However, even the value p m = 0.933 is larger than the 
upper limit of a possible liquid-solid tie line of p s = 0.912 [15] and p s = 0.904 
[13], respectively. The first result [15] was obtained by examining the behaviour 
of a local bond orientational order parameter as a function of the density. The 
second limit [13] was determined from a study of the pressure as a function 
of the density. Both simulations do not rule out a (one-stage or two-stage) 

2 Recent results concerning the liquid-solid structure of two-dimensional liquids can 
be found in Ref. [21]. 
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Fig. 1. FSS of the cumulant in the vicinity of p m . 

continuous transition, but give limits for a conventional first-order transition 
from the liquid to the solid phase. However, for such a first-order scenario the 
melting density p m (as determined from FSS of the positional order) must lie 
in the coexistence phase, i.e. below p s . Therefore, we can rule out a single 
first-order transition. 

From the scaling relation ip P o S 2 ~ L~ v we extract the critical exponent rj by a 
fit of all three system sizes. At p = 0.933 we get rj = 0.200(2), where the three 
points lie within statistical errors on a straight line. This result is incompatible 
with the prediction of the KTHNY theory of 1/4 < r)(p m ) < 1/3 [10]. Thus 
we can rule out the KTHNY scenario, because of the measured exponent rj. 
If our p m would be determined too low, this result remains unchanged. We 
performed also simulations in equilibrium at p = 1.0 (solid phase) to verify 
the scale invariance of U above p m . FSS yields rj = 0.0791(6). 

We now come to the dynamic relaxation of the hard disk system starting from 
the ordered state. Simulations of the short-time dynamics of systems with 
quasi- long-range order were performed for the 6-state clock model [22], the 
XY model [23,24] and the fully frustrated XY model [25,24]. From the scaling 
form of the second moment of the order parameter at or above p m 

V P os 2 (t,L) = b-^ p J(b- z t,b~ 1 L) (4) 



one obtains that the time dependence for sufficient large L is given by a power 
law of the form [5] 

^ P os 2 (t) ~ t-^ z . (5) 
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Accordingly, FSS analysis of the time-dependent cumulant 



U(t) = / p ° s W - 1 (6) 

(Vpos (*)) 

leads to 

U(t) ~ t d ' z . (7) 

Therefore, one can determine the dynamic critical exponent z from U (t) and 
then, with z in hand, the static exponent rj from the behaviour of ip p0 s 2 {t)- 111 
principle, one can also try to determine p m from the short-time dynamic, as 
it was done for a second order transition [26]. However, this determination is 
difficult in systems with a KT-like transition [23,5]. 

In the following we perform MC simulations at p m and in the solid phase 
(p = 1.0) to investigate the time evolution of il) vos 2 and U . For short-time sim- 
ulations we use the first definition of G (fixed direction), since large crystal 
tilting is not possible if we study only the early part of the evolution. How- 
ever, we checked that both definitions lead to similar results. Starting from 
the ordered state {jp pos = 1), i.e. with a perfect crystal with lattice spacing 
a, we release the system to evolve with the MC dynamic according to the 
Metropolis algorithm. We use systems of 8 2 , 16 2 , 32 2 , 64 2 and 128 2 hard disks 
and measure the observables up to 5000 MC sweeps. The average is taken over 
0(5000) samples for N = 128 2 and 0(300 000) samples for N = 8 2 . The criti- 
cal exponents are determined by least squares fits with a power law ansatz in 
the time interval t = (t m i n ,t max ). Statistical errors are calculated by dividing 
the data into different subsamples. Systematic errors are estimated by the re- 
sults of different system sizes and different time intervals, i.e. we examined the 
dependency of the slope from the fitted interval t = (t min , t max ) and number of 
particles N. The quoted error is a sum of the statistical and systematic error. 

In Fig. 2 (a) we plot the time evolution of ip vos 2 at the melting density for 
different system sizes in a double logarithmic scale. The figure shows that the 
power law behaviour starts after a microscopic time scale t m j C of approximately 
80 MC sweeps. For times up to 500 the difference between the systems with 64 2 
and 128 2 hard disks is negligible. Therefore, we omitted the data of iV = 128 2 
in the plot. Figure 2 (b) shows the time dependence of ^ pos 2 for N = 128 2 at 
both densities. In the time interval shown, the slope is nearly independent of 
time and yields r}/z = 0.0990(5) at p = 0.933 and rj/z = 0.0385(7) at p = 1.0, 
respectively. 

To determine z independently, we also measure the time evolution of the 
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Fig. 2. (a) Time evolution of the second moment of the positional order parameter 
V'pos 2 starting from the ordered state at p = 0.933. (b) as a function of time 

for p = 0.933 and p = 1.0. 

cumulant U. In Fig. 3 we show U(t) at the melting density for two different 
system sizes. From the slope we get z = 2.01(2). The analysis for p = 1.0 
yields z = 2.06(4). Thus we get rj = 0.199(3) at p = 0.933 and rj = 0.0794(29) 
at p = 1.0. These values coincide with the results from FSS within statistical 
errors, as can be seen from Table 1. Thus the data for rj are confirmed and 
new results for z are produced. 

In summary, we have performed short-time and equilibrium MC simulations 
of the hard disk model. The melting density p m and the critical exponent rj are 
determined from simulations in equilibrium. The results rule out the KTHNY 
scenario as well as a single first-order transition from the solid to the liquid 
phase. The short-time behaviour was used to extract the critical exponents r) 
and z at p = p m and in the solid phase. The values of 77 coincide with those 
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Fig. 3. Time dependent cumulant U(t) at p = 0.933 for N = 64 2 and N = 128 2 . 
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p Z T] 7] 



0.933 2.01(2) 0.199(3) 0.200(2) 

1.0 2.06(4) 0.0794(29) 0.0791(6) 
Table 1 

The critical exponents z and i] determined from the short-time dynamics of the 
system and the value of i] measured with FSS methods. 

from conventional FSS. Our simulations have shown that dynamic relaxation 
can also be used for this non-lattice model to measure the critical exponents. 
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